library(tidyverse)
library(data.table)
library(tidybulk)
library(ggpubr)

celline <-
fread('00_util_scripts/ref/cog-sanger-cell-line-rna-slim.csv.gz')

llps_ag <-
  'CXCR4
TNFRSF10B
EGFR
MET
PSMA
ERBB2
PDGFB
FZD1
FZD2
FZD5
FZD7
FZD8
LRP5
LRP6
ANPEP
CD33
CDC123
HAVCR2
CLEC12A
CD44
CD33
FTMT
TNFRSF14
TNFSF13B
TNFSF18
CD40LG
FGF1
FGF2
FGF4
FGF9
FGF10
TNFSF9
CD19
MS4A1' |> read_csv(col_names = 'gene')

ramos <- celline |>
  filter(model_name == 'Ramos' & symbol %in% llps_ag$gene)

ramos |>
  filter(fpkm > 0) |>
  mutate(symbol = fct_reorder(symbol, fpkm, .desc = TRUE)) |>
  ggplot(aes(symbol, fpkm, fill = symbol)) +
  geom_col() +
  geom_hline(yintercept = 20, linetype = 'dashed') +
  theme_pubr(x.text.angle = 45,legend = 'none') + labs_pubr() +
  labs(title = 'RNA expression in Ramos cell line',
       x = 'gene')

u937 <- celline |>
  filter(model_name == 'U-937' & symbol %in% llps_ag$gene)

u937 |>
  filter(fpkm > 0) |>
  mutate(symbol = fct_reorder(symbol, fpkm, .desc = TRUE)) |>
  ggplot(aes(symbol, fpkm, fill = symbol)) +
  geom_col() +
  geom_hline(yintercept = 20, linetype = 'dashed') +
  theme_pubr(x.text.angle = 45,legend = 'none') + labs_pubr() +
  labs(title = 'RNA expression in U-937 cell line',
       x = 'gene')

# process a20 dataset ------------
a20 <- read_delim('mission/GSE203162_GenewiseCounts.tsv.gz')

a20 <- a20 |>
  select(1:4) |>
  pivot_longer(3:4) |>
  mutate(GeneID = as.character(GeneID)) |>
  tidybulk(.sample = name, .transcript = GeneID, .abundance = value)

a20_scaled <- a20 |>
  identify_abundant() |>
  scale_abundance()

a20_scaled

## calc fpkm -----
a20 |> summarize(sum(value), .by = name)

a20_fpkm1 <- a20 |> filter(name == 'A201') |>
  mutate(fpkm = value * 1e9 / 49778949 / Length)

llps_ag_mm <- llps_ag$gene |>
  str_to_title() |>
  clusterProfiler::bitr(fromType = 'SYMBOL',
                        toType = 'ENTREZID',
                        OrgDb = 'org.Mm.eg.db')

a20_fpkm1 <- a20_fpkm1 |>
  rename(ENTREZID = GeneID) |>
  right_join(llps_ag_mm)

a20_fpkm1 |>
  filter(fpkm > 0) |>
  mutate(SYMBOL = fct_reorder(SYMBOL, fpkm, .desc = TRUE)) |>
  ggplot(aes(SYMBOL, fpkm, fill = SYMBOL)) +
  geom_col() +
  geom_hline(yintercept = 20, linetype = 'dashed') +
  theme_pubr(x.text.angle = 45,legend = 'none') + labs_pubr() +
  labs(title = 'RNA expression in A20 cell line',
       x = 'gene')

## calc fpkm on merged samples
a20_mrg <- a20 |>
  pivot_wider(names_from = 'name', values_from = 'value')

a20_mrg <- a20_mrg |>
  mutate(value = A201+A202, size = sum(value),
         fpkm = value * 1e9 / size / Length,
         ENTREZID = GeneID) |>
  select(ENTREZID, fpkm) |>
  right_join(llps_ag_mm)

a20_mrg |>
  filter(fpkm > 0) |>
  mutate(SYMBOL = fct_reorder(SYMBOL, fpkm, .desc = TRUE)) |>
  ggplot(aes(SYMBOL, fpkm, fill = SYMBOL)) +
  geom_col() +
  geom_hline(yintercept = 20, linetype = 'dashed') +
  theme_pubr(x.text.angle = 45,legend = 'none') + labs_pubr() +
  labs(title = 'RNA expression in A20 cell line',
       x = 'gene')

# sort GPCR expr in hpca -----------------
library(tidySummarizedExperiment)
hpca <- celldex::HumanPrimaryCellAtlasData()

gpcr_hs <- read_tsv('00_util_scripts/data/uniprot_human_gpcr.tsv.gz',
                    name_repair = 'universal')

hpca$label.main |>
  unique() |>
  str_subset('eratino')

hpca_kera <- hpca |>
  filter(label.fine == 'Keratinocytes') |>
  as_tibble()

hpca_kera_gpcr <- hpca_kera |>
  filter(.feature %in% gpcr_hs$Gene.Names..primary.)

hpca_kera_gpcr$.feature |> unique() |> length()

gpcr_hs$Gene.Names..primary. |> unique() |> length()

hpca_kera_gpcr |>
  mutate(.feature = fct_reorder(.feature, logcounts, .fun = mean)) |>
  slice_max(.feature, n = 45) |>
  ggplot(aes(.feature, logcounts)) +
  geom_boxplot() +
  theme_pubr(x.text.angle = 45) +
  expand_limits(y = 0) +
  labs(title = 'Top 10 GPCR mRNA expression in human keratinocytes',
       x = 'Gene', y = 'Normalized expression')

### GPCRs in 293T -----------
gpcrs <- query_uniprot_keyword('KW-0297')

hpa_cl <- fread('00_util_scripts/ref/rna_celline.tsv.gz', check.names = T)

gpcr_cl <- hpa_cl |>
  filter(Gene.name %in% gpcrs$symbol)

gpcr_mean <- gpcr_cl |>
  summarise(avg_nTPM = mean(nTPM), .by = Gene.name)

gpcr_mean

gpcr_cl |>
  filter(Cell.line == 'HEK293', nTPM > 0) |>
  left_join(gpcr_mean) |>
  filter(nTPM > avg_nTPM) |>
  as_tibble() |>
  mutate(rank = row_number(-nTPM)) |>
  slice_max(nTPM, n = 20) |>
  ggplot(aes(rank, nTPM)) +
  geom_point() +
  geom_label_repel(aes(label = Gene.name), nudge_y = 10) +
  geom_area(aes(y = avg_nTPM), alpha = .5) +
  annotate(geom = 'text', x = 4, y = 2, color = 'white',
           label = 'Average nTPM in all cell lines') +
  theme_bw() +
  scale_x_reverse() +
  labs(x = 'Rank', title = 'Top 20 expressed GPCRs in HEK293 cell line',
       subtitle = 'Human protein atlas')

## lyko2020 ---------
lyko <- read_rds('00_util_scripts/data/GSE130973_seurat_analysis_lyko.rds')
source('00_util_scripts/mod_seurat.R')
library(harmony)
library(SingleR)

lyko |>
  ggplot(aes(x = 'x', percent.mito)) + geom_violin()

lyko <- lyko |>
  UpdateSeuratObject()

lyko |>
  VlnPlot('ADRB2')

Idents(lyko) |> table()

DimPlot(lyko)

lyko$integrated_snn_res.0.4 |> table()

lyko <- lyko |>
  mutate(seurat_clusters = integrated_snn_res.0.4) |>
  mark_cell_type_singler(hpca, new_label = 'hpca_main')

lyko_kera <- lyko |>
  filter(str_detect(hpca_main, 'Kerat'))

lyko_kera |>
  get_abundance_sc_long(features = gpcr_hs$Gene.Names..primary.) |>
  select(.feature, .abundance_RNA) |>
  mutate(gene = fct_reorder(.feature, .abundance_RNA, ExpMean)) |>
  slice_max(gene, n = 13990) |>
  ggplot(aes(gene, .abundance_RNA)) +
  geom_violin() +
  stat_summary(geom = 'crossbar', fun = 'ExpMean', color = 'red', width = .5) +
  theme_pubr(x.text.angle = 45) +
  labs(title = 'Top 10 GPCR mRNA expression in human keratinocytes',
       subtitle = 'GSE130973 (scRNA-seq)',
       x = 'Gene', y = 'Normalized expression')

# mouse kera -----------
library(Seurat)
library(tidyseurat)
sobj_kera <- read_rds('mission/FPP/zww_sa_mice/data/aureus_mice_skin_kera.rds')

gpcr_mm <- gpcr_hs$Gene.Names..primary. |>
  str_to_title()

srt_avrg <- sobj_kera |>
  AddMetaData('one', col.name = 'all') |>
  AverageExpression(group.by = 'all', features = gpcr_mm)

srt_avrg |>
  as.data.frame() |>
  as_tibble(rownames = 'gene') |>
  mutate(gene = fct_reorder(gene, RNA)) |>
  slice_max(RNA, n = 10) |>
  ggplot(aes(gene, RNA)) +
  geom_col() +
  theme_pubr() +
  labs(title = 'Top 10 GPCR mRNA expression in mouse keratinocytes',
       x = 'Gene', y = 'Normalized expression')

top10mm <- srt_avrg |>
  as.data.frame() |>
  as_tibble(rownames = 'gene') |>
  mutate(order = fct_reorder(gene, RNA)) |>
  slice_max(RNA, n = 10)

sobj_kera |>
  get_abundance_sc_long(features = c(top10mm$gene, 'Adra1a')) |>
  mutate(order = fct_reorder(.feature, .abundance_RNA, ExpMean)) |>
  ggplot(aes(order, .abundance_RNA)) +
  geom_violin() +
  stat_summary(geom = 'crossbar', fun = 'ExpMean', color = 'red', width = .5) +
  theme_pubr() +
  labs(title = 'Top 10 GPCR mRNA expression in mouse keratinocytes',
       x = 'Gene', y = 'Normalized expression')

## GSE246720 ----
zhou2023 <- read_delim('00_util_scripts/data/GSE246720_L14_RNA_seq.txt.gz')

zhou_ordered <- zhou2023 |> select(2,3,6) |>
  filter(GeneSymbol %in% gpcr_mm) |>
  reframe(w1 = ExpMean(WT1), w2 = ExpMean(WT2), .by = GeneSymbol) |>
  pivot_longer(2:3) |>
  mutate(gene = fct_reorder(GeneSymbol, value))
  
top10zhou <- zhou_ordered |>
  slice_max(gene, n = 20) |>
  pull(GeneSymbol) |> unique()

zhou_ordered |>
  filter(gene %in% c('Adra1a','Adrb2',top10zhou)) |>
  ggplot(aes(gene, value)) +
  geom_boxplot() +
  theme_pubr() +
  labs(title = 'Top 10 GPCR mRNA expression in mouse keratinocytes',
       subtitle = 'GSE246720 keratinocyte bulk-seq',
       x = 'Gene', y = 'Normalized expression')


hpa.sc.cls <- read_delim('hpa/rna_single_cell_type_tissue.tsv')

hpa.sc.cls$Tissue |> unique()

pbmc.ln.ccr6 <- hpa.sc.cls |>
  filter(`Gene name` == 'CCR6' & Tissue %in% c('pbmc','lymph node'))

pbmc.ln.ccr6 |>
  filter(`Cell type` != 'mixed cell types') |>
  ggplot(aes(`Cell type`, nTPM, color = `Cell type`)) +
  geom_boxplot() +
  facet_grid(~Tissue, scales = 'free_x', space = 'free_x') +
  labs(title = 'CCR6 expression in human protein atlas') +
  theme_pubr(x.text.angle = 45)
